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In this note, we assess the predictive validity of stochastic dynamic causal modeling 
(sDCM) of functional magnetic resonance imaging (fMRI) data, in terms of its 
ability to explain changes in the frequency spectrum of concurrently acquired 
electroencephalography (EEG) signal. We first revisit the heuristic model proposed in 
Kilner et al. (2005), which suggests that fMRI activation is associated with a frequency 
modulation of the EEG signal (rather than an amplitude modulation within frequency 
bands). We propose a quantitative derivation of the underlying idea, based upon a neural 
field formulation of cortical activity. In brief, dense lateral connections induce a separation 
of time scales, whereby fast (and high spatial frequency) modes are enslaved by slow 
(low spatial frequency) modes. This slaving effect is such that the frequency spectrum of 
fast modes (which dominate EEG signals) is controlled by the amplitude of slow modes 
(which dominate fMRI signals). We then use conjoint empirical EEG-fMRI data — acquired 
in epilepsy patients — to demonstrate the electrophysiological underpinning of neural 
fluctuations inferred from sDCM for fMRI. 



Keywords: dynamic causal modeling, neural noise, EEG, fMRI, effective connectivity, neural field, separation of 
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INTRODUCTION 

In addition to localizing brain regions that encode specific sen- 
sory, motor, or cognitive processes, contemporary neuroimaging 
research tries to understand how information is transmitted 
within brain networks (Sporns, 2010). Ultimately, the ambition 
is to understand the nature of the information that neuronal 
populations pass on to each other. This speaks to the notion of 
functional integration, which views cognitive function as result- 
ing from information exchange within brain networks (Tononi 
et al., 1994). This means one has to understand the directed influ- 
ence or effective connectivity that brain systems (e.g., cortical areas, 
neuronal populations, or single neurons) exert on each other. 
Analyzing effective connectivity rest on models that formalize 
assumptions about how neuronal systems are wired and how they 
respond to different stimuli. These models are then used to inter- 
pret brain responses measured using, e.g., functional magnetic 
resonance imaging (fMRI) or magneto-/electroencephalography 
(M/EEG). We refer to Valdes-Sosa et al. (2011) for a compre- 
hensive review on effective connectivity, and its relationship with 
influence and causality. 

In the context of M/EEG, detailed biophysical knowledge 
about the activity of neural ensembles has led the community to 
develop and validate dynamical models that can describe macro- 
scale dynamics in great detail (cf. Deco et al, 2008; Coombes, 
2010; Bressloff, 2012). Most of these models are inspired by 



approaches in statistical physics based on the notion of a mean 
field, i.e., the idea that interactions among ensembles of neu- 
rons can be captured by summary statistics (i.e., moments of 
the relevant distribution). Dynamic causal modeling [DCM, see 
Daunizeau et al. (2011) for a recent review] for M/EEG embeds 
these models into a formal (Bayesian) statistical framework that 
allows for parameter estimation and model comparison when 
analyzing evoked or induced responses. In this context, the 
approach has proven successful in exploiting the realism of these 
biophysical models to capture experimental effects of cognitive 
manipulations in terms of network plasticity (e.g., Garrido et al., 
2007 or Moran et al, 201 1). 

Even more established is the use of DCM for fMRI data, for 
which the approach was originally proposed by Friston et al. 
(2003), building on previous advances in the biophysics of hemo- 
dynamic processes (Buxton et al., 1998). Interestingly, fMRI 
paradigms elicit neural responses that unfold on a much slower 
time scale (seconds) than that typical of electrophysiological mea- 
surements (milliseconds). In contrast to hemodynamic processes, 
relatively little is known about the biophysical underpinning of 
these neuronal responses, which has precluded the use of detailed 
(realistic) models of neural dynamics in DCM for fMRI. Instead, 
the slow (and widespread) dynamics of interacting brain regions 
are captured by phenomenological models that are much simpler 
than those employed for M/EEG data. Even though these models 
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serve as a reference when evaluating fMRI effective connectivity 
methods (Smith et al., 2011), their simplicity may impose lim- 
itations on the ensuing inferences (but see, for example, David 
et al., 2008 or Brodersen et al, 2011). This is why we have 
previously proposed variational Bayesian approaches that can 
deal with some of the (unavoidably simplifying) model assump- 
tions by introducing random fluctuations of physiological states 
(Friston et al, 2008, 2010; Daunizeau et al, 2009b). Initial appli- 
cations of this "stochastic" version of DCM [stochastic dynamic 
causal modeling (sDCM)] have established some face validity in 
the context of fMRI data analysis. For example, Friston et al. 

(2011) showed how sDCM can be used in combination with 
post-hoc model comparison (Friston and Penny, 2011) to explore 
large model spaces, and Li et al. (2011) examined the smooth- 
ness of the underlying neural fluctuations. In Daunizeau et al. 

(2012) , we provide a comprehensive comparison of determinis- 
tic and stochastic DCM, in terms of parameter estimation and 
model comparison. In particular, we showed how accounting for 
random effects on the system's dynamics can improve network 
identification by exploiting the decorrelation of neural time series 
induced by neural noise. 

However, sDCM for fMRI has its own methodological chal- 
lenges and still lacks solid experimental validation. On the one 
hand, the ability of sDCM for fMRI to discriminate between 
the respective contributions of neuronal noise and measurement 
noise in the observed fMRI signal depends on modeling assump- 
tions (Roebroeck et al., 2011; Daunizeau et al., 2012). On the 
other hand, it is yet unclear how to formally relate the neu- 
ral states of DCM for fMRI to electrophysiological measures of 
activity. In this note, we use EEC data to provide empirical evi- 
dence for the predictive validity of sDCM for fMRI data. We 
therefore place this work in the context of the ongoing debate 
regarding the hemodynamic correlate of the EEG [see Rosa et al. 
(2010a) for a recent critical review of this literature]. Here, we 
appeal to neural field theory (Wilson and Cowan, 1972, 1973; 
Lopes da Silva et al, 1973; Nunez, 1974; Amari, 1975, 1977; 
Freeman, 1975; Jirsa and Haken, 1996) to revisit the heuristic 
hemodynamic correlate of EEG signals originally proposed in 
Kilner et al. (2005). In essence, this hypothesis suggests that the 
fMRI signal is correlated with the frequency modulation of the 
EEG (where larger fMRI signals reflect a shift of the EEG spec- 
trum toward higher frequencies). This is based on the idea that 
slow dynamics underpinning fMRI responses reflect the instan- 
taneous frequency of oscillating modes of activity (Cabral et al., 
2011; Friston et al., 2011). Here, we show how slow macroscopic 
modes of activity that emerge from intrinsic (within region) con- 
nections may modulate the frequency spectrum of fast dynam- 
ical modes that contribute to the EEG signal. Guided by this 
notion, we obtain estimates of neural dynamics from fMRI data 
using sDCM and ask whether these estimates predict the spec- 
tral behavior of concurrent EEG data — in terms of the temporal 
changes in its center frequency (above and beyond physiological 
confounds). This enables us to establish the predictive valid- 
ity of sDCM for fMRI data, in relation to electrophysiological 
responses. 

This paper comprises three sections. In the first, we present 
the neural field model that motivates a separation of time scales 



and the ensuing heuristic hemodynamic correlate of the EEG. 
The second section is an empirical demonstration that serves 
to evaluate the predictive validity of stochastic DCM for fMRI 
for the ensuing electrophysiological correlates. We close with 
a discussion of the implications of our findings in the final 
section. 

MODEL: REVISITING THE HEURISTIC HEMODYNAMIC 
CORRELATE OF EEG 

In this section, we first introduce a neural field model of a cor- 
tical region, and show how its spatiotemporal pattern of activity 
can be decomposed into canonical dynamic modes that have dis- 
tinct time scales. In brief, slow macroscopic neural states control 
the rate of change of fast components. This motivates both the 
separation of time scales that underlies the generative model of 
sDCM for fMRI (Friston et al., 2011) and the heuristic hemody- 
namic correlate of EEG (Kilner et al, 2005). In brief, this section 
serves as a partial theoretical derivation of sDCM for fMRI, from 
the perspective of biologically-informed (neural field) models of 
macroscopic electrophysiological activity. 

NEURAL FIELDS AND SEPARATION OF TIME SCALES 

We start with a description of the dynamics of single neurons 
within an ensemble, which we describe in terms of the proba- 
bility density function p(q) of neuronal post-synaptic membrane 
potentials (PSP) q. We assume that neurons fire an action poten- 
tial if their PSP surpasses a depolarization threshold x- This 
means that (H(q — x)) represents the average (mean field) firing 
rate over the ensemble, where H is the Heaviside step function 
and ( ) is the expectation under the ensemble density. Using the 
central limit theorem, one can show that the average firing rate 
over the ensemble is determined mainly by the first two moments 
of the PSP ensemble density p(q) (Liley et al., 2002; Marreiros 
etal.,2008): 

(H( S - X )> = f°°p(s)dq « S(T)) = — J- - 

J x l + exp(-p(i) - x)) 

n = <s> (D 

P~ 2 = 4«£-ri) 2 ) 

where S is a sigmoid function and p is proportional to the 
inverse of the standard deviation (dispersion) of depolarization 
within the ensemble. This is a statistical (mesoscopic) summary 
of neuronal ensemble activity, corresponding roughly to a corti- 
cal macrocolumn (about 100,000 neurons and 1 mm 2 of cortex). 
For our modeling purposes, this summary statistics is now taken 
to the macroscale, where each ensemble or column has a posi- 
tion r on the 2D-cortical manifold. At this macroscopic scale, 
states like depolarization can be regarded as a continuum or field 
(Jirsa and Haken, 1996; Liley et al., 2002; Coombes et al, 2007; 
Daunizeau et al, 2009a; Pinotsis et al, 2012), which is a function 
of space and time: v\(t) — > r\(r, f). This means we can describe 
the local average input and output activity of neurons using 
the depolarization field r\(r, f) and the sigmoid function. The 
spatiotemporal dynamics of the neural field essentially depends 
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on how local ensembles influence each other — as described 
below. 

We now consider a region of the cortical manifold that is 
defined by its continuous domain D, over which the neural field 
r\(r, t) is defined. The field is considered to result from the 
temporal convolution of pre-synaptic input u,(r, t) as follows: 



/oo 
j(t - tVo, t)dt' 
-oo 



: / (g) (x(r, f) 



(2) 



7(f) 



j Y ^exp(-^) f>0 
0 t < 0 



where / is the synaptic impulse response function, x is a time- 
constant, and y controls the maximum PSP following an action 
potential. The maximum PSP and time-constant summarize the 
microscopic properties of the ensemble (Jansen and Rit, 1995). 
Because action potentials propagate with a finite conduction 
velocity, the mean rate of arrival of pre-synaptic impulses can 
be expressed as a time-delayed integral (spatio temporal convo- 
lution) of local firing rates: 



,x(r, t)= f f G(|r - r'|, t - f')S(ti(r / , t'))dr'dt' 

J-oo JD 



(3) 



G<g)Sor|(r, f), 



where D is the domain, over which the neural field is defined, and 
G is a homogeneous Green function that acts as a spatiotemporal 
convolution operator. Let r = | r — r' | be the distance between two 
points with positions r and r'. Below, we follow Bojak and Liley 



(2010) and parameterize the Green function as follows: 

H{t), (4) 



G° 1 , / r 2 + v 2 t 2 

G(r, t) = -r 1 exp 

4jt a 2 F V 2avt 



where the Heaviside step function endows Equation (3) with 
temporal causality and G° is the total number of synaptic con- 
nections. This dispersive propagator assumes a Gaussian falloff in 
the density of synaptic connections and accounts for a certain 
variability in velocities (v is the mode of the velocity distribution 
at "large distance," i.e., when r/a — » oo). The spatial scale of lat- 
eral connectivity is controlled by the spatial dispersion a. Figure 1 
depicts the dependence of the dispersive propagator with respect 
to space and time (see also Liley et al., 2012 — Figure 3). 

In addition to its plausible form, this homogeneous Green 
function has proven very useful, in that it has an analytical Fourier 
transform. Applying the Fourier transform back and forth to 
Equations (2) and (3) yields the following partial differential 
equation (see Bojak and Liley, 2010): 



v 9 
2a + dt 



IX (r, t) 



GV 
2a 



S(tl(r, 0) 



(5) 



X*—, + 2x— + 1 ii Or, t) = x Y |x(r, t). 
at at 



The first partial differential equation describes the diffusion of 
action potentials of spikes (through the Laplacian operator) and 
derives from the lateral connectivity kernel in the Green function. 
The second ordinary differential equation derives from Equation 
(2) and describes synaptic transmission; i.e., how propagated pre- 
synaptic firing is accumulated locally to produce depolarization. 



16 



log-density of connections 




30 40 50 
time (msec) 
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FIGURE 1 | Neural field model: the dispersive propagator. This figure 
depicts the dispersive propagator of Equation (4) as a function of both time 
(x-axis) and distance (y-axis), in terms of the density 2nrG(r, f) of synaptic 
connections that are reached by an action potential emitted at r= t = 0 
(on a logarithmic scale). Note that the 2-nr scaling arises because of radial 



symmetry in 2D (Bo]ak and Liley, 2010). The white line shows the average 
distance, at which activity is propagated as a function of time. The black 
lines are contour lines of the connection density. These can be used to 
eyeball how fast the density attenuates (e.g., along the average distance 
path). 
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The solution to Equation (5) can be decomposed into canon- 
ical dynamical modes, by projecting the neural field on the set 
of eigenvectors (r) of the Laplacian operator defined over the 
domain D (see Daunizeau et al., 2009a). These spatial modes solve 
the Laplacian eigenvalue problem, i.e.,: V 2 w^'(r) = ^-k w ( r )' 
where the distribution of eigenvalues X is such that: Xjt < Xo = 0. 
They form a complete and orthonormal basis function set for the 
neural field, allowing for an eigen-decomposition of the field, as 
follows: 



r|(r, t) = J2 z ? W w{k) W 

k 

Li(r, t) = J2zf ) (t)w (k Ht), 



(6) 



through the coupling function c ( z i° This slaving effect is such 
that the frequency profile of fast modes (k > 0) is controlled by 
the amplitude of slow modes (k = 0). 

MODELING LOCAL ACTIVATION IN TERMS OF EEG FREQUENCY 
MODULATION 

Equation (8) expresses the Jacobian of eigenmode dynamics as a 
function of the fundamental mode of activity zj 0 ' within a given 
brain region, over which the neural field is defined. This means 
that slow fluctuations in the fundamental mode zj should be 
associated with changes in the EEG frequency spectrum P(co), 
which derives from the Laplace transform of Equation (8) (see 
Appendix 2): 



here, the dynamics of the spatial modes are defined as: 



w W (r)r|(r, t)dr 



(7) 



* (l) (r)|i(r, t)dr. 



Furthermore, it can be shown that the dynamics of the spatial 
modes (Equation 5) can be approximated as third-order ordinary 
differential equations that are coupled through the fundamental 
mode k = 0 (see Appendix 1): 



■Jk) 



z 2 



z w + 0(z 2 ) 



(8) 



(,<»), !£,(,») (,- S (,»>)) 



2a 

here, Z3 = 'z\, and c (zf^ is a non-linear function of the fun- 
damental mode Zj°^ that mediates the coupling between modes 
describing firing rate and depolarization. Roughly speaking, each 
mode decays at a rate that is proportional to a\>X k . As a con- 
sequence, high propagation velocities in large cortical regions 
dissipate the spatial modes quickly, with the exception of the fun- 
damental mode, which has a zero eigenvalue. More generally, 
smooth eigenmodes (low spatial scales) will be associated with 
slow dynamics (see, e.g., Schultze-Kraft et al, 2011). Note that 
the natural time scale of the fundamental mode can be very slow, 
if the conduction velocity v becomes small compared to the dis- 
persion of lateral connections a (i.e., if v/a — »• 0). In this case, 
the fundamental mode z®' will decay so slowly that it will pre- 
dominate over other stable (higher spatial frequency) modes z^\ 
which disappear as soon as they are created. This means that 
Equation (8) effectively instantiates a separation of time scales, 
where the slowest mode "enslaves" the fast modes (Haken, 1983) 



P(w) oc 



( to +£(1-0^)) 
(»+! + !) + ? 



y - 



(9) 



where k c is the cut-off order induced by the skull's low electrical 
conductivity. Here, the summation is over the fast modes (up to 
k c = 16), and we have assumed that EEG responses reflect post- 
synaptic potentials, i.e., z\ ^ . Similarly to Kilner et al. (2005), let 
us now define the EEG center frequency w as the first moment of 
the (normalized) frequency power spectrum P(w): 



P(o>) 



P((0) 



/ 



P(w)c2w 



(10) 



/ 0L)P(a))cZa) 



By construction, normalized frequency power P(u)) is the 
proportion of total power attributed to any given frequency of 
interest. Thus, one can think of it in terms of a probability den- 
sity function over frequencies, whose first-order moment would 
be the center frequency. Intuitively, an increase in the center 
frequency d> signals a shift of frequency power toward higher fre- 
quencies. Note that the center frequency is not the main peak in 
the frequency spectrum. Instead, variations of a> over time quan- 
tify frequency modulations of the EEG signal. The key idea here 
is that the EEG center frequency w is a function of the slow eigen- 
mode dynamics zf\ and its modulation over time thus follows 
the slow eigenmode dynamics zf \ i.e., d> = w ( z i°')- 

Figure 2 depicts the frequency modulation of the EEG signal 
induced by changes in the slow eigenmode z f\ The frequency 
power spectrum P(w) given in Equation (9) was evaluated using 
the parameter values given in Table 1 below. One can see that 
when the slow eigenmode tends toward the action-potential firing 
threshold (z| 0 ' — >■ x)> there is both a power increase in the low 
frequencies and a decrease in the high frequencies (cf. increase 
in the delta/theta/alpha band, and decrease in the beta/gamma 
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FIGURE 2 | Neural field model: EEG frequency modulation. This 
figure shows the effect of changes in the mean membrane 
depolarization fundamental eigenmode z{ 0) on the frequency content of 
the EEG signal. Upper-left: EEG frequency spectrum P(<o) (z-axis) as a 
function of frequency io (x-axis) and eigenmode amplitude zj 0 ' (y-axis). 
Upper-right: Normalized frequency spectrum P(w) (cf. Equation 10, 



same format). Lower-left: Magnitude of standard EEG frequency bands 
(y-axis; blue: delta, green: theta, red: alpha, magenta: beta, violet: 
gamma) as a function of the eigenmode amplitude zf* (x-axis). 
Lower-right: Centre frequency u> (y-axis) as a function of the 
eigenmode amplitude zj 0 ' (x-axis). Note that here, the action potential 
firing threshold x is 30 mV. 



band). This shift toward lower frequencies can be seen most 
clearly on the normalized frequency power P(oo). As a conse- 
quence, the center frequency 6j increases as the slow eigenmode 
gets further away from the action potential firing threshold. In 
our example, this induces a change of about 10 Hz in the EEG 
center frequency. Note that the precise numerical value of the 
center frequency a> (as well as its susceptibility to the eigenmode 
zj 0 ^) depends upon the neural field's parameters (cf. Table 1). 
In fact, using more realistic neural fields models (accounting 
for, e.g., excitatory and inhibitory interneurons, see Pinotsis 
et al, 2012) would profoundly change the shape of the fre- 
quency power spectrum. However, the qualitative effect of neural 



activation upon the center frequency is qualitatively invariant 
to such changes. We will comment on this in the "Discussion" 
section. 

THE LINK BETWEEN EEG TIME-FREQUENCY RESPONSE AND fMRI 
TEMPORAL RESPONSE 

In this paper, we assume that the slow macroscopic dynamics x'"' 
underlying fMRI data in a given region of interest follow the fun- 
damental eigenmode zj 0 ' of its membrane depolarization cortical 
field. The motivations for this assumption include the following: 

• Its relative simplicity. 
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Table 1 | Neural field parameters. 



Constant 


Physical meaning 


Value 


Y 


Maximum post-synaptic impulse response 


8mV 


T 


Post-synaptic impulse response decay rate 


4 ms 


G° 


Total number of synaptic connections 


2000 


a 


Spatial decay rate of synaptic connections 


3 mm 


P 


Inverse of the within-ensemble PSP 


0.54 mV" 1 




standard deviation 




X 


Action potential firing threshold (w.r.t. PSP 


30 mV 




steady-state) 




V 


Average conduction velocity 


3m/s 



following Langevin equation for stochastic DCM (see Daunizeau 
et al, 2012): 



These parameter values were taken from Pinotsis et al. (2012). 



• Sensitivity to input versus output spiking activity: there is 
a consensus now on the fact that mean membrane depolar- 
ization (as opposed to e.g., firing rate) drives hemodynamic 
changes [see Logothetis and Wandell (2004) for a review on 
the relationship between local field potentials and the BOLD 
signal). 

• Energy budget: both the slow time scale and the low spatial 
frequency of the fundamental eigenmode make it likely to be 
the main driver of induced BOLD changes, when compared 
with faster (and higher spatial frequency) eigenmodes. This is 
because blood flow increases in local vascular trees are induced 
by coherent spatio-temporal summation of neural activity (see 
e.g., Daunizeau et al., 2010; Rosa et al., 2010a,b; Schmuel, 
2010). 

• It relates to theoretical accounts on critical slowing of high- 
dimensional dynamical systems, which has recently been used 
to motivate the rate of change of neural states in stochastic 
DCM for fMRI (Friston et al, 2010). 

We will further comment on this point in the "Discussion" 
section. 

In the previous section, we have described how a separation 
of time scales could arise from local interactions at the level of 
a cortical region. It turns out that this result can be generalized 
to a set of coupled local neural fields, in that the fundamen- 
tal mode mediates the extrinsic (between-region) coupling. In 
other words, the influence different neural fields exert on each 
other expresses itself through the coupling of their respective 
fundamental modes, which locally drive higher order modes (as 
in Equation 8). Deriving the precise dynamical properties of 
coupled neural fields is beyond the scope of the present study. 
However, one can invoke the slaving effect when considering cou- 
pled (distal) brain regions, the network properties of which shape 
the temporal dynamics of local fundamental eigenmodes zf\ 
This means that long-range connections control the frequency 
profile of brain regions responses (see e.g., Bojak et al, 2011), 
through their effect on their respective fundamental eigenmodes 
Zj 0 ' . Since we did not consider coupled neural fields in the above 
treatment, we replace their (unknown) evolution function with 
a first-order Taylor expansion on slow neural states, yielding the 



■An) 



Ax (n) + v 



(11) 



where 8 W = {A,-,} are neuronal parameters that measure extrinsic 
coupling strengths. Here, we have neglected non-linear interac- 
tion terms that act as gating factors on the network connections 
x^ (cf. Stephan et al, 2008). The influence of the fast modes 
on the motion of slow modes is expressed via fluctuations v, 
which are comparatively much faster. This allows us, in the con- 
text of sDCM for fMRI, to treat v as stochastic (neural) noise and 
place priors p(\\m) on its generalized motion (see below). 

In contradistinction, typical evoked EEC responses are driven 
by a mixture of slow and fast modes. Note that the relative con- 
tribution of the slow and fast modes is imbalanced (toward fast 
modes), simply because the summation span all eigenmodes up to 
cut-off order k c . In fact, typical evoked EEC responses disappear 
within a second or so of peristimulus time; i.e., their limit fre- 
quency has an order of magnitude similar to the fMRI sampling 
rate (about 1 Hz). 

Taken together, we expect to see concurrent changes in the slow 
neural states x^ 1 ' that drive BOLD and in the EEC power spec- 
trum P(w). However, from the above section, we know that the 
underlying relationship is non-linear and depends upon intrin- 
sic properties of local neural fields, such as the size of active 
brain regions and the spatial decay rate of synaptic connections 
(cf. Equation 9). In addition, the respective contribution of each 
region to the scalp EEG power spectrum is virtually unknown, 
as this would require solving the so-called EEG inverse prob- 
lem (see Rosa et al, 2010a,b). These issues conspire with the 
unavoidable simplifications of our model to make the predic- 
tion of the full EEG frequency spectrum practically irrelevant. 
In the following, we will thus focus on the EEG center fre- 
quency d), which captures global changes in the frequency spec- 
trum. For the sake of simplicity, in the remainder of this work, 
we will invoke a first-order Taylor expansion of the following 
form: 



(*«) 



w(0) + 



3x« 



x (n) +0(x 2 ) 



(12) 



where the gradient 9a>/9x (n) | Q measures the susceptibility of the 
center frequency to the slow neural states. This quantity will be 
directly estimated from experimental data (along with the inter- 
cept d»(0), which acts as a confound). Testing for Equation (12) 
effectively bypasses the dependency of the EEG power spectrum 
to the respective contribution each local neural field included 
in the network, as well as their biophysical properties, and only 
retains the qualitative prediction of the model. 

Summary 

In this section, we have used the neural field formulation of local 
neural activity to demonstrate an emergent separation of tempo- 
ral scales, in which the slow fluctuations of modes or patterns of 
depolarization enslave faster modes. Crucially, the amplitude of 
slow modes (which we associate with fMRI signals) and modulate 
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the frequencies of fast modes (which we associate with EEG sig- 
nals). This dynamical behavior is based upon a plausible mean 
field approximation to macroscopic neuronal activity on the 
cortical surface and, within the setting, provides a formal ver- 
ification of the heuristic in Kilner et al. (2005), linking fMRI 
signals to changes in the center frequency of electrophysiologi- 
cal signals. In the next section, we will use stochastic DCM to 
obtain estimates of macroscopic (slow) neural fluctuations from 
fMRI data. We then test whether these estimates can predict fluc- 
tuations in concurrently measured EEG data — in terms of its 
center frequency, as suggested by the neural field treatment above 
(Equations 9-12). 

EMPIRICAL DATA: RESTING STATE AND GENERALIZED 
SPIKE AND WAVE (GSW) ACTIVITY 

In this section, we conduct an sDCM analysis of data from 
patients with epilepsy who exhibited generalized spike and wave 
(GSW) activity, while undergoing concurrent EEG/fMRI mea- 
surements in the scanner. We first describe the EEG/fMRI acqui- 
sition and conventional activation analysis. We then describe the 
sDCM analysis of fMRI data and ensuing validation using EEG 
data. 

EEG/fMRI ACQUISITION AND ACTIVATION ANALYSIS 

These data were part of a previous neuroimaging study of idio- 
pathic generalized epilepsy (Hamandi et al, 2006; Vaudano et al., 
2009), whose acquisition protocol and data pre-processing are 
briefly summarized here. 

Ten-channel EEG (10-20 system) was recorded using MR- 
compatible equipment, along with bipolar electrocardiogram. 
After filtering and artifact correction (Krakow et al, 2000), 
the onsets and offsets of GSW activity were identified by two 
experts (see Vaudano et al., 2009 for details). Seven hun- 
dred T2* -weighted single-shot gradient echo echo-planar images 
(TE = 40 ms, TR = 3 s, 21 interleaved axial slices of 5 mm 
thickness, FOV = 24 x 24 cm 2 , 64 x 64 matrix) were acquired 
over a 35min session with a 1.5 Tesla MRI scanner (Horizon 
EchoSpeed, General Electric). Patients were asked to rest with 
their eyes closed and to keep still. FMRI data were pre-processed 
using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/). EPI time series 
were realigned and spatially smoothed with an 8 mm FWHM 
isotropic Gaussian kernel and spatially normalized to the stan- 
dard anatomical space. For each patient, a general linear model 
(GLM) was constructed to test for the presence of regional GSW- 
related BOLD changes. Periods of GSW activity were modeled 
as blocks, beginning at GSW onset and terminating at their off- 
set. The GSW regressors were then convolved with the canonical 
hemodynamic response function (plus temporal and dispersion 
derivatives) before inclusion in the GLM. In addition, we included 
slow drifts (Fourier basis function set), as well as motion-related 
effects (head and eye movements), cardiac-related effects (see 
Liston et al., 2006) and scan-nulling regressors (modeling inter- 
scan motion events larger than 0.2 mm, cf. Lemieux et al., 2007) 
as confounding factors. 

On average, the GLM design matrices contained about 100 
confounds regressors. Significant positive and negative GSW- 
related BOLD responses were identified by means of an F-contrast 



on the GSW regressors for the nine patients included in this 
study. SPMs were thresholded atp < 0.05 (FWE whole-brain cor- 
rected) to define three regions of interest (ROIs), which were 
implicated in the initiation and termination of GSW discharges 
in all patients: thalamus, prefrontal cortex (PFC), and precuneus. 
A summary time series was derived for each ROI by computing 
the first eigenvariate of all suprathreshold voxel time series within 
a 10 mm of the ROI centers. The time series were adjusted for all 
confounding effects included in the GLM analysis. 

VALIDATING THE sDCM MACROSCOPIC NEURAL STATES ESTIMATES 
WITH EEG: METHODS 

In addition to the neural evolution function given in Equation 
(11), DCM for fMRI requires the specification of an additional set 
of hemodynamic states that couple neural dynamics to observed 
BOLD signal changes: 

xW-Mf-K/^-l) 

x« = To(^ 2 ~* ' ) e 3 (13) 

e x 2 '4 [l-d-EoY 2 1 

J_ V L - gU-a)4 /a 

to E 0 

- V /- 

here, x- n > represents regional neural activity, whose dynamics 
is given in Equation (11). Equation (13) expresses changes in 
hemodynamic states x^ 1 ', as a response to a neural perturba- 
tion x^ n) and as a function of hemodynamic parameters 6™ = 
{k s , Kf, a, to. Eq, eo} (for details, see the appendix in Stephan 
et al, 2008). Finally, one has to specify the observation mapping 
from hemodynamic states x 1 -' 1 ^ to observed local BOLD changes y 
(Stephan et al., 2007): 

y = V 0 (4.3 v 0 E 0 TE (l - e*< J + e 0 r 0 E 0 TE (l - e*4 -*> J 
+ (l-e 0 )(l-^" ) ))+8, (14) 

where e is an additive measurement noise, and the parameters are 
defined in Table 2 below. 

In brief, the predicted data y depend non-linearly on the 
unknown model variables © = {x, 0}, through Equations (11), 
(13), and (14). This model (as well as statistical assumptions 
about measurement noise e) is encoded in the likelihood func- 
tion p (y 1 0, ni). Priors p (0 | m ) specify our assumptions about 
the magnitude of state noise, evolution parameters and obser- 
vation parameters — where the prior means of the hemodynamic 
parameters of Equations (13) and (14) are given in Table 2: 

Inverting the generative model m means (1) approximat- 
ing the conditional density p (0|y, m) of unknown variables 0 
given the set of sampled measurements y and (2) quantifying 
the model evidence p (y | m ) . Non-linearities in the generative 
model eschew exact analytical solutions to this problem, which 
is finessed using variational approaches that rest on optimizing a 
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Table 2 | Prior means of hemodynamic parameters and 
acquisition-related constants for fMRI at 1.5 T. 



Constant 


Physical meaning 


Value 




Vacnni atnrw cinna riorai; rata 
VdbUUIIdlLHy biyild UcUdy [die 


0 65 Hz 




\/acnni atnru cinnu Toonnaw rata 
VdbUUIIdlUly blylld ccUUdLK Idle 


0 41 Hz 


to 


Mean transit time 


2s 


« 


Vessel stiffness 


0.32 


Eo 


Oxygen extraction fraction at rest 


0.34 


Vo 


Venous volume fraction at rest 


4 


»0 


Frequency offset 


40.3 Hz 


TE 


EPI echo time 


0.04 s 


r 0 


Intravascular relaxation rate 


25 Hz 


so 


Ratio of intra- and extra-vascular signals 


1 



Note f/iaf 9 < ' 1) = (k s , K f , a, to, Eo, eo} are estimated from the data (i.e., have 
a positive prior variance), whereas the remaining parameters are treated as 
constants. 

free-energy lower bound F(q) to the model evidence, with respect 
to an approximate conditional density q(&): 



F(q) = (Inp (0 \m) + hip (y |0, m) - hiq (0)) 
= \np (y \m) - D KL (<?(©); p(0 |y, m)) , 



(15) 



where Dkl is the Kullback-Leibler divergence and the expecta- 
tion ( ) q is taken under q. From Equation (15), maximizing the 
functional F(q) with respect to q minimizes the Kullback-Leibler 
divergence between q(@) and the exact posterior p(0 |y, m). 
This decomposition is complete in the sense that if q(&) = 
p(@ |y, m), then F(q) = lnp(y \m). Typically, the iterative max- 
imization of free energy proceeds under the Laplace approxi- 
mation, where the approximate posterior q(&) ~ p(0 |y, m) is 
assumed to have a Gaussian form (see Friston et al, 2007). We 
refer to Daunizeau et al. (2009b) for details about the application 
of the variational Bayesian approach to stochastic DCM. 

Following model inversion as described above, inference on 
hidden states is based on the conditional density q(x) and thus 
depends upon the generative model m. Usually, in DCM, infer- 
ence on states or parameters is preceded by model selection. 
However, our focus was on the validity of sDCM for fMRI 
with regard to predicting the spectral behavior of concurrently 
acquired EEG data — based on the idea that slow macroscopic 
neural states x'"' control the EEG frequency modulation. We 
therefore endowed the network with full reciprocal connectivity 
(no zero entry in matrix A of Equation 1 1 ) and simply applied the 
above variational Bayesian approach to derive the first moment 
x = (x'"' ) of the conditional density q(x) . In the following, x thus 
refers to the sDCM estimate of the slow neural states, given fMRI 
time series y (and under model m). 

To assess the relationship between EEG frequency modulation 
and neural states as estimated by sDCM, we constructed a GLM 
Hi, whose dependent variable was the trajectory w of the EEG 
center frequency sampled at each fMRI time sample and whose 
independent variables were the neural states estimated by sDCM 



plus confounds (cf. Equation 12): 



Hi : co = [x Xo] 



PoJ 



+ e 



(16) 



where x are the sDCM conditional estimates, Xo contains the con- 
founds (see below), P and P 0 are unknown regression coefficients, 
and e are i.i.d. residuals (with zero mean and unknown variance). 
Parameters P effectively capture the susceptibility of the EEG cen- 
ter frequency with respect to the slow neural states (cf. gradient 
9w/9x^"'| 0 in Equation 12). They lump together the contribu- 
tions of the multiple regions that constitute the network as well 
as biophysical properties of the local neural fields. Parameters P 0 
weight the contribution of the confounding factors to changes 
in the EEG center frequency. We then tested whether the sDCM 
neural states predicted EEG frequency modulation above and 
beyond what could be explained using the confounds. To this 
end, we used Bayesian model comparison at the group level 
(Stephan et al., 2009) to quantify the evidence in favor of the full 
model Hi , relative to a reduced model Ho that only included the 
confounds: 



H 0 : w = X 0 Po + e 



(17) 



The assessment of sDCM predictive validity thus depends on 
the definition of the confounds, i.e., uncontrolled sources of vari- 
ations in the EEG center frequency. In order to investigate the 
robustness of the statistical relationship between the slow eigen- 
mode dynamics x and the EEG center frequency trajectory <o, we 
have defined two different sets of confounds: 

• Slow drifts: this assumes that the EEG center frequency under- 
goes slow and unspecific variations due to e.g., skin conduc- 
tance changes, residual eye blinks artifacts, etc Under 

this assumption, we set up Xo similarly to the minimal con- 
founds of a typical fMRI activation analysis, i.e., a Fourier basis 
function set (16 first harmonics; see Figure 6). 

• Full confounds: this generalizes the intuition that confounds 
on the EEG center frequency are similar to those included in 
the fMRI activation analysis. For example, cardiac and respi- 
ratory effects, as well as body movements, could be modeled 
directly, as potential sources of (nuisance) variations. Under 
this assumption, we set up Xo identically to the full con- 
founding matrix used for each subject fMRI activation analysis 
(cf. section "EEG/fMRI Acquisition and Activation Analysis" 
above). 

In brief, testing for the contribution of sDCM neural estimates 
x on variations of the center frequency w, above and beyond the 
full confounds is the most conservative of the two alternatives. 
However, it will be the least sensitive of both models as well; recall 
that P has only mroi entries, where mroi = 3 is the number of 
ROIs included in the sDCM analysis (to be compared with the 
dimension of ft 0 ). In other words, these two sets of confounds 
can be motivated on principled grounds. We will thus perform a 
statistical comparison between Ho and Hi that marginalizes over 
the two sets of confounds. 
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VALIDATING NEURAL STATES ESTIMATE OBTAINED BY sDCM 
AGAINST EEG 

Figure 3 shows the SPM of significant GSW-related BOLD 
changes (above and beyond physiological confounds), for the first 
patient of the study (F-test, p < 0.05, whole-brain corrected). 
This was used to define the three ROIs included in the sDCM 
analysis (Figure 4). 

One can see that most of the variance of the BOLD response is 
explained by the stochastic DCM. Focusing on peri-GSW periods, 
one can see that the conditional density q(x) shows high tempo- 
ral variability, with increased uncertainty during head motion. 
This is expected, as scan-nulling confound regressors effectively 
prevent sDCM from deriving any information from fMRI data 
about hidden states during these periods. It is interesting to note 
that the network connectivity coefficients (elements of the A 
matrix in Equation 11) are small (about 0.02 Hz). This means 
that the average coupling of these regions (over the entire record- 
ing) is rather weak. Importantly, with the (partial) exception 
of the inhibitory intrinsic self-connections, the posterior cor- 
relation matrix (Figure 4) shows that most system parameters 
are identifiable and not confounded by hemodynamic parame- 
ter estimation. Generally, Figure 4 summarizes the motivation for 
this work — is the temporal variability of sDCM estimates of neu- 
ronal activity an artifact of measurement noise, or does it have an 
electrophysiological underpinning? 

The cross-validation of neural state estimates by sDCM 
from fMRI data against EEG data requires the extraction of 
the time-dependent center frequency co for each channel of 
each subject, according to Equation (10). This is depicted in 



Figure 5, for the EEG channel C4 of the same patient in 
Figures 3, 4. 

Recall that comparisons between the experimental and theo- 
retical profiles of the frequency spectrum show differences that 
prevent meaningful quantitative fitting of Equation (9). For 
example, the model cannot account for peaks in the alpha band 
that appear in occipital channels (at least for the neural field 
parameter values of Table 1; results not shown). However, it is 
reassuring to see that the variations of the EEG frequency pro- 
file over time appear as small perturbations around a largely 
invariant average power spectrum. More precisely, the appear- 
ance of GSW activity does not dramatically distort its shape, 
but rather induces frequency modulations of small amplitude. 
This is important, since this affords face validity to the quali- 
tative prediction of our neural field model. Interestingly, there 
is a clear change in the normalized frequency power begin- 
ning at the onset of the second GSW crisis. This leads to a 
marked decrease in the EEG center frequency (the first GSW event 
seems to be associated with a similar but weaker effect). We will 
briefly comment on GSW-related frequency modulations in the 
discussion. 

Figure 6 summarizes the relationship between slow macro- 
scopic sDCM states and the EEG frequency modulation of chan- 
nel C4 for the same patients as in previous figures. Even though 
one can see that the center frequency d> is very well fitted by 
the model (upper row of Figure 6: coefficient of determination: 
R 2 = 0.59), one might think that most of the frequency modu- 
lation could be explained by confounds. For example, eyeballing 
the trajectory of the observed EEG center frequency reveals clear 
low frequency tends that are likely to be captured by slow drifts 
confounds. To disclose the specific contribution of neural state 
estimates, we have plotted the data and model prediction after 
removing the confounds (adjusted data). One can see that, when 
considering slow drifts confounds, sDCM neuronal estimates 
explain EEG frequency modulations of about 10 Hz. This is a 
similar order of magnitude than center frequency variations cap- 
tured by slow drifts (cf. parameter estimates on lower-right panel 
of Figure 6). Note that the corresponding correlation coefficient 
is about 0.49 (adjusted R 2 = 0.24), which indicates a reasonably 
good fit (above and beyond confounds). 

Figure 7 summarizes the result of Bayesian model compar- 
ison, in terms of the log-Bayes factor LBF = logp(w \H\ ) — 
logp((o |Ho), across all EEG channels for the same patient as 
before. One can see that all EEG channels show strong evi- 
dence in favor of the full model (Hi), whereas there is evidence 
against Hi only for EOG and ECG channels. This suggests that 
the neuronal fluctuations estimated by stochastic DCM predict 
the EEG frequency modulation beyond slow drifts confounds. 
For completeness, we also show the results of classical inference 
(F-test with significance level p = 0.05, full model Hi against 
the null Ho). As this test was repeated independently for each of 
the 30 channels, we applied a Bonferroni correction for multi- 
ple comparisons across channels. One can see that the profile of 
F-statistics and log-Bayes factors (LBFs) across channels is very 
similar. 

Recall that the results summarized in Figures 6, 7 are only 
shown to illustrate the analysis. In fact, it is important to note 




FIGURE 3 | Absence seizure analysis: regions of interest. This figure 
summarizes the standard SPM activation results of a single epileptic 
patient with (petit mal) absence seizure (the same subject as in Figure 1). 
Significant (whole brain FWE-corrected, p < 0.05) positive and negative 
GSW-related BOLD responses were identified by means of an F-contrast 
on the GSW regressors. The color bar indicates the range of displayed 
F values. 
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FIGURE 4 | Absence seizure analysis: sDCM analysis. This figure 
summarizes the sDCM analysis of a single subject fMRI data (same 
patient as in Figure 1). Upper-left: Time series of the observed and 
fitted BOLD signal in prefrontal cortex (PFC) as a function of time. The 
red shaded area indicates peri-GSW activity. Upper-right: Observed 
(y-axis) versus fitted (x-axis) BOLD signal in PFC. Lower-left: Macroscopic 
sDCM neural states in the three regions of interest (PF, prefrontal cortex; 



pC, precuneus; Th, thalamus) during peri-GSW activity, as a function of 
time. The red bar indicates strong head motion, which was modeled 
using scan-nulling regressors (see main text). Lower-right: sDCM 
conditional density on model parameters (left: first-order moment and 
ensuing network connectivity, right: posterior correlation matrix). 
Note: identifiability issues between pairs of parameters manifest as high 
posterior correlations. 



that they (1) depend upon the definition of the confounds (here 
slow drifts) and (2) do not capture inter-individual variability. 
The question of whether the accuracy with which the neural states 
inferred by sDCM are able to predict the frequency modulations 
of the EEG signal (above and beyond confounds) generalizes at 
the group level is addressed below. 

This Bayesian model comparison for the single patient shown 
in Figure 7 was then repeated for all subjects, for both sets of con- 
founds (namely: slow drifts and full confounds). First, we pooled 
evidence over EEG channels for each subjects. This was simply 
done by summing the log- model evidences over EEG channels 
(within-subjects fixed effect. The ensuing distribution of LBFs 
across subjects — shown in Table 3 below — indicates that most 
subjects show very strong evidence in favor of Hi, irrespective of 
the particular set of confounds. 

To confirm this, we ran a random-effect analysis (Stephan 
et al, 2009) to assess the prevalence of Hi and Ho at the group 
level. Note that the two sets of confounds ("slow drifts" and "full 
confounds") induce the four following models: 



• H} ormsj : Hi with slow drifts 

• h[ ' : Hi with full confounds 

• HQ drlfts) : Ho with slow drifts 

• Hg full ' ) : Ho with full confounds. 

The comparison of interest is in terms of whether or not adding 
the sDCM slow neural states adds anything to the prediction of 
the EEG center frequency, above and beyond confounds. Thus, 
we applied group-level family inference (Penny et al, 2010) to 
evaluate the evidence in favor or against Hi , irrespective of the 
definition of the confounds. We have thus partitioned the set of 
four models above into two families Hi and Ho that gather both 
sorts of confounds, as follows: Hi = j j/j drlfts ) | H{ full) J and Ho = 

| Hg drifts) , Hq 11 ^ J . Deriving the relative evidence of Hi against Ho 

effectively marginalizes over the different confounds. First of all, 
we checked that there was a random effect at all, i.e.,: a difference 
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FIGURE 5 | Absence seizure analysis: derivation of the EEG frequency 
modulation. This figure depicts the extraction of the EEG center frequency in 
a single subject (same patient as in Figure 1). Left: The EEG set-up used during 
the recording session is shown superimposed on the brain and skin surfaces 



(sensor C4 is highlighted). Right: Normalized square root power of windowed 
Fourier transform (z-axis) of the EEG traces of sensor C4 (x-axis: scanning 
time, y-axis: instantaneous frequency). The blue line shows the center 
frequency <u as a function of scanning time f (cf. Equation 10 in the main text). 



between the respective prevalence of all models within the pop- 
ulation. We found that the posterior probability of all models 
being equally frequent was p = 0.03. This Bayesian omnibus 
test indicates that the distribution of Bayes factors shown in 
Table 3 is unlikely to be due to chance. Figure 8 summarizes the 
group-level family inference results. One can see that the esti- 
mated prevalence of family Hi is about 0.82. In addition, the 
exceedance probability of Hi (versus Ho) was 0.991. In other 
words, we can be 99% confident that sDCM estimates of neu- 
ronal fluctuations predict concurrent EEG frequency modulation 
above and beyond physiological confounds (in more than half the 
population). 

For completeness, we have also performed bayesian 
group-level model comparison, conditional on the definition 
of confounds. The results match those of the family-inference 
analysis above. In brief, the exceedance probability of £fj drlfts) 
(versus HQ dnfts ') was 0.997, and the exceedance probability of 
jj(Ml) ( versus ff tf"") ) was 0.990. However, the bayesian omnibus 
test for these analyses was not as conclusive (slow drifts: p = 0. 1 1, 
full confounds: p = 0.14). Taken together, this means that there 
is statistical evidence in favor of Hi (versus Ho), irrespective of 
the definition of the confounds. This concludes our quantitative 
assessment of the predictive validity of sDCM for fMRI with 
regard to concurrently obtained EEG measures. 

DISCUSSION 

Introducing random neural fluctuations in DCM for fMRI was 
originally motivated by the need to estimate hidden states in 
the absence of experimental perturbations. In this note, we have 
assessed the predictive validity of sDCM estimates of slow macro- 
scopic neural states with respect to EEG frequency modulation. 
We have revisited the heuristic hemodynamic correlate of the EEG 
proposed in Kilner et al. (2005), using a neural field formulation 
to show how slow modes of local electrophysiological activity 



(which we assume drive changes in the fMRI signal) enslave 
the frequency spectrum of fast modes that contribute to the 
EEG signal. From this theoretical perspective, we tested whether 
slow neural estimates obtained by sDCM could explain concur- 
rent EEG frequency modulation above and beyond physiological 
confounds. 

To our knowledge, this work is the first attempt to pro- 
vide empirical evidence that neural fluctuations inferred using 
stochastic DCM from fMRI time series have an electrophysi- 
ological underpinning. From an experimental perspective, this 
endorses the use of sDCM for studying spontaneous fluctuations 
that are beyond experimental control. Examples of this include — 
but are not limited to — multistable perception (Hesselmann et al, 
2010) and epileptic events (e.g., absence seizures, cf. Vaudano 
et al, 2009). In this context, the added-value of sDCM is that it 
allows for comparing quantitative hypotheses about the propa- 
gation and/or the emergence of these phenomena within brain 
networks. 

The acute reader might question the relevance of the modeling 
part of our paper, in the sense that quantitative predictions of the 
model were not used in the subsequent data analysis. However, 
we believe this model may contribute to the debate regarding 
the hemodynamic correlate of the EEG. First, it bears a number 
of interesting differences, with respect to the heuristic model of 
Kilner et al. (2005): 

• Here, the main prediction, namely that fMRI activation is 
associated with a frequency modulation of the EEG signal, 
only relies on the assumption that BOLD changes are mostly 
driven by the fundamental (slow and smooth) eigenmode of 
the membrane potential field. This differential sensitivity of 
EEG and fMRI with respect to higher-order canonical eigen- 
modes bypasses most modeling assumptions of the heuristic 
model. 
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FIGURE 6 | Absence seizure analysis: EEG results (1). This figure 
summarizes the results of the analysis on the EEG frequency modulation 
in a single subject (same patient as in Figure 5) at channel C4. Upper-left: 

Observed (blue dashed line) and predicted (black plain line) frequency 
modulation (y-axis) of the best EEG channel (C4) as a function of scanning 
time (x-axis). For this patient, the two GSW episodes are indicated by the 
orange arrows. Upper-middle: Observed (y-axis) versus predicted (x-axis) 



frequency modulation for EEG channel C4. Upper-right: Matrix of confounds 
Xo included in the analysis (slow drifts). Lower-left: Observed (blue dashed 
line) and predicted (black plain line) frequency modulation (y-axis) of the best 
EEG channel (C4) as a function of scanning time (x-axis), after adjustment for 
confounds. Lower-middle: Adjusted (y-axis) versus predicted (x-axis) 
frequency modulation for EEG channel C4. Lower-right: GLM parameter 
estimates, ± one standard deviation (orange: p, gray: p 0 ). 



• Here, the relationship between activation and the EEG 
center frequency is not monotonic. More precisely, the 
center frequency decreases as one approaches the action 
potential threshold, either from hyperpolarized states 
(activation) of from depolarized states (de-activation). In 
brief, this form of critical slowing is qualitatively equiv- 
alent to the heuristic model, above the action potential 
threshold. 

Second, we believe that the above qualitative theoretical predic- 
tion is very robust to model assumptions. First, we checked that 
the effect persisted in the context of our model, when changing 
the parameters that control the biophysical properties of the neu- 
ral field (results not shown). Second, the effect is a byproduct 
of local propagation of activity over the field. In brief, smooth 
local patterns emerge in densely connected brain regions, the — 
slow — dynamics of which are shaped by effective connectivity 
between remote regions. Each local eigenmode is then expected 
to control the acceleration or slowing down of the response of its 
respective brain region. In contradistinction to e.g., quantitative 
predictions about frequency power spectra, we anticipate 



this qualitative prediction to hold, irrespective of models' 
specificities. 

Let us now consider the (modest) theoretical contribution 
of this work in relation to the separation of time scales. We 
believe that — if properly extended — established biophysical mod- 
els of fast electrophysiological responses may be a good starting 
point for understanding slow dynamical modes at macroscopic 
spatial scales that, presumably, drive fMRI responses. The neu- 
ral field treatment we have developed in this note is a first 
step toward a mechanistic understanding of the macroscopic 
dynamics that emerge at the time scale of seconds or minutes. 
However, this model is limited in many ways. First, we did not 
account for extrinsic connectivity between regions (coupled neu- 
ral fields). This is known to strongly impact on the stability of 
spatiotemporal brain dynamics (Knock et al., 2009), and will 
thus be critical for predicting the dynamical repertoire of slow 
macroscopic modes of activity (cf. point above). For example, 
Bojak et al. (2011) show that increasing the long-range cou- 
pling between remote regions strongly modulates the frequency 
response of the target region, in terms of the relative power in 
low- and high- frequency bands. Second, we did not separate the 
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FIGURE 7 | Absence seizure analysis: EEG results (2). This figure 
summarizes the results of the analysis on the EEG frequency modulation 
in a single subject (same patient as in Figures 5, 6) across channels. 
Upper panel: Value of the F-statistic (y-axis) when testing for the 
significance of the contribution of the sDCM neural states dynamics in 
the EEG frequency modulation, as a function of EEG channels (x-axis). 
Red stars indicate channels that pass the 5% false positive rate 



threshold (with Bonferroni correction of the multiple comparisons across 
channels). Upper-right inset: Associated p-value across EEG channels 
(the red line indicates the corrected 5% threshold). Lower panel: Log 

Bayes factor log p(H-\ \y) — logp(Ho \y) showing evidence in favor of Hi 
versus Ho across EEG channels (see main text). Red lines indicate 
posterior probabilities of model Hi of 95% (upper line) and 5% (lower 
line), respectively. 



Table 3 | Log-Bayes factors (LBF) of H-\ versus Ho, for all subjects, and 
both sets of confounds. 



Subject 


LBF (Hj against 


LBF (Hj against 




slow drifts) 


full confounds) 


1 


670.9 


198.6 


2 


204.3 


103.9 


3 


51.4 


-4.8 


4 


92.9 


23.3 


5 


17.9 


185.4 


6 


4.4 


35.3 


7 


-1.8 


65.0 


8 


1.2 


28.6 


9 


101 .4 


37.2 



A positive LBF indicates evidence in favor of H\ (jLBFI> 3: significant effect 
at the subject-level, Kass and Raftery, 1995). Highlighted cells depict the two 
subjects that show relative evidence in favor of Ho- 



respective contributions of excitatory and inhibitory subpopu- 
lations (Pinotsis et al, 2012). Of particular importance here is 
the fact that inhibito-excitatory networks can possess stable limit 
cycles (Seung et al., 1995). This means that local neural fields 
can behave as oscillators or resonators. This relates to the work 
of Cabral et al. (2011) who examined the role of local network 
oscillations in the emergence of slow temporal synchrony between 
remote brain areas. The critical thing here is that the time it takes 
for one resonator to influence another resonator can be much 
slower than the local resonance frequencies. In addition, owing to 
the autonomous nature of local dynamics, the effect of between- 
regions coupling might be better understood in terms of changes 
in the resonance frequency of the target region. Taken together, 
these considerations suggest that challenging theoretical work still 
needs to be undertaken to establish the dynamical repertoire of 
coupled slow macroscopic modes of brain activity. 

We take the fact that the sDCM neural states contribution 
was significant for all EEG channels (above and beyond slow 
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FIGURE 8 | Absence seizure analysis: EEG results (3). This figure 
summarizes the results of the analysis on the EEG frequency modulation at 
the group level (across subjects). Left: Posterior estimates of the expected 
frequency of families Ho and within the population (± one posterior standard 



deviation). Right: Posterior density (y-axis) over the expected frequency 
(x-axis) of family within the population. The red shaded area indicates the 
mass of probability that is beyond 50% prevalence, yielding an exceedance 
probability of about 99%. 



drifts) as a compelling validation of the conditional estimates 
of neuronal activity from sDCM. Importantly, this also provides 
some empirical evidence in favor of the fMRI to EEG mapping 
itself. This supports the results of Rosa et al. (2010b), which is, 
to our knowledge, the only experimental study that tested pre- 
dictions of the heuristic model of Kilner et al. (2005). In this 
study, the authors experimentally controlled the frequency con- 
tent of the EEG signal, and showed that the BOLD response in 
a group subjects was better explained by frequency modulation 
than by amplitude modulation of the EEG. More specifically, Rosa 
et al. (2010b) showed that to correctly model BOLD variations 
it was important to use a normalized frequency spectrum, as in 
Equation (10). Here, we approach this issue from the other side 
and use fMRI to predict the EEG frequency spectrum. Note that 
this is only an indirect validation of the neurovascular coupling 
model, as we did not compare it with other quantitative scenar- 
ios. Among these is the idea that BOLD changes correlate with 
amplitude modulations within frequency bands of interest in the 
EEG spectrum. Typically, the BOLD signal is considered a fil- 
tered version of the EEG alpha power (see e.g., De Munck et al., 
2008). More generally, this idea has been applied to each and every 
EEG frequency band with relative success [see Laufs et al. (2008) 
and references therein]. As an illustration, Lachaux et al. (2007) 
have found that EEG gamma band modulations co-localize with 
BOLD variations. However, although power fluctuations of dif- 
ferent EEG frequency bands are known to be mutually highly 
correlated, multi frequency models are usually not used in such 
studies (De Munck et al., 2009). In addition to compromising the 
specificity of the results, this prevents a direct comparison with 
our and related work (cf. Rosa et al., 2010b). 

Another limitation of our data analysis is that we neglected 
the unavoidable non-linearity implicit in the relationship w = 

w (zj 0 ^ (cf- Equations 9-12 and Figure 2) when testing for the 
contribution of slow neural states to modulations of the EEG cen- 
ter frequency. This means we may have not taken full advantage of 
our neural field formulation. However, this does not detract from 



our quantitative assessment of the predictive validity of sDCM for 
fMRI data. 

When examining the GSW-related EEG frequency modula- 
tions, we noted that the frequency spectrum did not change 
markedly during GSW activity. This might seem surprising, if one 
thinks of GSW as sudden bursts of activity within the network, 
associated with clearly recognizable 3 Hz transient waves of EEG 
activity (Destexhe and Sejnowski, 2001). Closer inspection of the 
EEG frequency profile confirmed the appearance of a peak at 3 Hz 
during GSW-activity (not shown). However, many other (faster) 
processes actually contribute to the overall frequency spectrum, 
which typically drags the center frequency toward 20 Hz. This is 
the main reason why we had to include as many confounding fac- 
tors as possible, when assessing the specific contribution of the 
sDCM predictors. Importantly, physiological confounds (cardiac 
activity, head movements, etc.) span most of the slow frequencies 
in the EEG signal, which potentially masks interesting slow neu- 
ral processes. This means that slow frequency EEG signals with 
a neural origin are confounded with physiological noise. Note 
that there are notorious exceptions to this detectability issue, for 
example readiness potentials (Jahanshahi and Hallett, 2003) that 
are thought to reflect the preparation of motor activity and have a 
typical time scale of a few seconds. Nevertheless, our results indi- 
cate that it may be possible to detect the contribution of slow 
neural processes to the scalp EEG signal, through their indirect 
effect on the frequency modulation of fast modes. 

In conclusion, we have provided empirical evidence support- 
ing the predictive validity of sDCM for fMRI data. In future work, 
we will examine the empirical and theoretical relationships that 
exist between slow and fast dynamical consequences of effective 
connectivity as can be accessed via fMRI and EEG, respectively. 

ACKNOWLEDGMENTS 

The authors would like to thank Dr. Dimitris Pinotsis for use- 
ful discussions. K. J. Friston was funded by the Wellcome Trust. J. 
Daunizeau was funded by an ERC grant. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



January 2013 | Volume 6 | Article 103 | 14 



Daunizeau et al. 



Predictive validity of sDCM 



REFERENCES 

Amari, S.-I. (1975). Homogeneous 
nets of neuron-like elements. Biol. 
Cybern. 17,211-220. 

Amari, S.-I. (1977). Dynamics of pat- 
tern formation in lateral-inhibition 
type neural fields. Biol. Cybern. 27, 
77-87. 

Bojak, I., and Liley, D. T. (2010). 
Axonal velocity distributions in 
neural field equations. PLoS Comp. 
Biol. 6:el000653. doi: 10.1371/jour- 
nal.pcbi. 1000653 

Bojak, I., Oostendorf, T. R, Reid, A. 
T., and Kotter, R. (2011). Towards 
a model-based integration of 
co-registered electroencephalogra- 
phy/functional magnetic resonance 
imaging data with realistic neu- 
ral population meshes. Philos. 
Transact. A Math. Phys. Eng. Sci. 
369, 3785-3801. 

Bressloff, P. C. (2012). Spatiotemporal 
dynamics of continuum neural 
fields. /. Phys. A Math. Theor. 
45:033001. doi: 10.1088/1751-8113/ 
45/3/033001 

Brodersen, K. H., Schofield, T. 
M., Leff, A. P., Ong, C. S., 
Lomakina, E. I., Buhmann, J. 
M., et al. (2011). Generative 
embedding for model-based clas- 
sification of fMRI data . PLoS 
Comput. Biol. 7:el002079. doi: 
10.1371/journal.pcbi.l002079 

Buxton, R. B., Wong, E. C, and 
Franck, L. R. (1998). Dynamics 
of blood flow and oxygenation 
changes during brain activation: the 
Balloon model. Magn. Res. Med. 39, 
855-864. 

Cabral, J., Hugues, E., Sporns, O., 
and Deco, G. (2011). Role of local 
network oscillations in resting-state 
functional connectivity. Neuroimage 
57, 130-139. 

Coombes, S. (2010). Large-scale neu- 
ral dynamics: simple and complex. 
Neuroimage 52, 731-739. 

Coombes, S., Venkov, N. A., Shiau, 
L., Bojak, L, Liley, D. T. J., and 
Laing, C. R. (2007). Modeling 
electrocortical activity through 
local approximations of integral 
neural field equations. Phys. Rev. E 
Stat. Nonlin. Soft Matter Phys. 76, 
1539-1547. 

Daunizeau, J., David, O., and Stephan, 
K. E. (2011). Dynamic Causal 
Modelling: a critical review of 
the biophysical and statistical 
foundations. Neuroimage 58, 
312-322. 

Daunizeau, J., Kiebel, S. J., and Friston, 
K. J. (2009a). Dynamic causal mod- 
eling of distributed electromagnetic 
responses. Neuroimage 47, 590-601. 

Daunizeau, J., Kiebel, S. J., and Friston, 
K. J. (2009b). Variational Bayesian 



identification and prediction of 
stochastic dynamic causal models. 
Physica D 289, 2089-2118. 

Daunizeau, J., Laufs, H., and Friston, 
K. J. (2010). "EEG-fMRI informa- 
tion fusion: biophysics and data 
analysis," in EEG-fMRI- Physiology, 
Technique and Applications, eds C. 
Mulert and L. Lemieux {Berlin, 
Heidelberg: Springer), 511-526. 
ISBN: 978-3-540-87918-3. 

Daunizeau, J., Stephan, K. E., and 
Friston, K. J. (2012). Stochastic 
dynamic causal modelling of fMRI 
data: should we care about neural 
noise? Neuroimage 62, 464-481. 

David, O., Guillemain, I., Baillet, S., 
Reyt, S., Deransart, C, Segebarth, 
C, et al. (2008). Identifying neu- 
ral drivers with functional MRI: 
an electrophysiological validation. 
PLoS Biol 6:e315. doi: 10.1371/jour- 
nal.pbio.0060315 

Deco, G., Jirsa, V. K., Robinson, P., 
Breakspear, M., and Friston, K. 
(2008). The dynamic brain: from 
spiking neurons to neural masses 
and cortical fields. PLoS Comp. 
Biol. 4:el000092. doi: 10.1371/jour- 
nal.pcbi.l 000092 

De Munck, J. C. (1988). The potential 
distribution in a layered spheroidal 
volume conductor. /. Appl. Phys. 64, 
464^70. 

De Munck, J. C, Goncalves, S. I., Faes, 
T. J., Kuijer, J. P., Pouwels, P. J., 
Heethaar, R. M., et al. (2008). A 
study of the brain's resting state 
based on alpha band power, heart 
rate and fMRI. Neuroimage 42, 
112-121. 

De Munck, J. C, Goncalves, S. I., 
Mammoliti, R., Heethaar, R. M., 
and Lopes da Silva, F. H. (2009). 
Interactions between different 
EEG frequency bands and their 
effect on alpha-fMRI correlations. 
Neuroimage 47, 69-76. 

Destexhe, A., and Sejnowski, T. J. 
(2001). Thalamocortical Assemblies. 
Oxford, UK: Oxford University 
Press. 

Freeman, W. J. (1975). Mass Action 
in the Nervous System: Examination 
of the Neurophysiological Basis of 
Adaptive Behavior through the EEG, 
1st Edn. New York, NY: Academic 
Press. 

Friston, K. J., Harrison, L., and 
Penny, W. D. (2003). Dynamic 
causal modelling. Neuroimage 19, 
1273-1302. 

Friston, K. J., Li, B., Daunizeau, J., and 
Stephan, K. E. (2011). Network dis- 
covery with DCM. Neuroimage 56, 
1202-1221. 

Friston, K. J., Mattout, J., Trujillo- 
Barreto, N., Ashburner, J., and 
Peeny, W. (2007). Variational free 



energy and the Laplace approxima- 
tion. Neuroimage 34, 220-234. 

Friston, K. J., and Penny, W. (2011). 
Post-hoc Bayesian model compari- 
son. Neuroimage 56, 2089-2099. 

Friston, K. J., Stephan, K. E., Li, 
B., and Daunizeau, J. (2010). 
Generalised filtering. Math. 
Probl. Eng. 2010:621670. doi: 
10.1155/2010/621670 

Friston, K. J., Trujillo-Barreto, N. J., 
and Daunizeau, J. (2008). DEM: a 
variational treatment of dynamical 
systems. Neuroimage 42, 849-885. 

Garrido, M. I., Kilner, J., Kiebel, S. J., 
and Friston, K. J. (2007). Evoked 
brain responses are generated by 
feedback loops. Proc. Natl. Acad. Sci. 
U.S.A. 104, 20961-20966. 

Haken, H. (1983). Synergetics, an 
Introduction: Nonequilibrium Phase 
Transitions and Self-Organization 
in Physics, Chemistry, and Biology, 
3rd Rev. enl. ed. New York, NY: 
Springer- Verlag. 

Hamandi, K., Salek-Haddadi, A., Laufs, 
H„, Liston, A., Friston, K., Fish, 
D. R., et al. (2006). EEG-fMRI 
of idiopathic and secondarily gen- 
eralized epilepsy. Neuroimage 31, 
1700-1710. 

Hesselmann, G., Sadaghiani, S., 
Friston, K. J., and Kleinschmidt, 
A. (2010). Predictive coding or 
evidence accumulation? False 
inference and neuronal fluctu- 
ations. PLoS ONE 5:e9926. doi: 
10.1371/journal.pone.0009926 

Jahanshahi, M., and Hallett, M. (eds.). 
(2003). The Bereitschaftspotential, 
Movement-Related Cortical 
Potentials. Springer Series. New York, 
NY: Kluwer Academic / Plenum 
Publishers. ISBN: 0306474077. 

Jansen, B. H., and Rit, V. G. (1995). 
Electroencephalogram and visual 
evoked potential generation in a 
mathematical model of coupled cor- 
tical columns. Biol. Cybern. 73, 
357-366. 

Jirsa, V, and Haken, H. (1996). Field 
theory of electromagnetic brain 
activity. Phys. Rev. Lett. 77, 960-963. 

Kass, R. E., and Raftery, A. E. (1995). 
Bayes factors. /. Am. Stat. Assoc. 90, 
773-795. 

Kilner, J. M., Mattout, J., Henson, 
R., and Friston, K. J. (2005). 
Hemodynamic correlates of 
EEG: a heuristic. Neuroimage 28, 
280-286. 

Knock, S. A., Mcintosh, A. R., Sporns, 
O., Kotter, R., Hagmann, P., and 
Jirsa, V. K. (2009). The effects 
of physiologically plausible con- 
nectivity structure on local and 
global dynamics in large scale brain 
models. /. Neurosci. Methods 183, 
86-94. 



Krakow, K., Allen, P. J., Lemieux, L., and 
Fish, D. R. (2000). Methodology: 
EEG-correlated fMRI. Adv. Neurol. 
83, 187-201. 

Lachaux, J. P., Fonlupt, P., Kahane, P., 
Minotti, L., Ho mann, D., Bertrand, 
O., et al. (2007). Relationship 
between task-related gamma oscil- 
lations and BOLD signal: new 
insights from combined fMRI 
and intracranial EEG. Hum. Brain 
Mapp. 28, 1368-1375. 

Laufs, H., Daunizeau, J., Carmichael, 
D. W., and Kleinschmidt, A. (2008). 
Recent advances in recording 
electrophysiological data simulta- 
neously with magnetic resonance 
imaging. Neuroimage 40, 515-528. 

Lemieux, L., Salek-Haddadi, A., Lund, 
T. E., Laufs, H., and Carmichael, 
D. (2007). Modelling large motion 
events in fMRI studies of patients 
with epilepsy. Magn. Reson. Imaging 
25, 894-901. 

Li, B., Daunizeau, J., Stephan, K. E., 
Penny, W., and Friston, K. (2011). 
Stochastic DCM and generalized fil- 
tering. Neuroimage 58, 442-457. 

Liley, D. T. J., Cadusch, P. J., and Dafilis, 
M. P. (2002). A spatially continuous 
mean field theory of electrocortical 
activity. Network 13, 67-113. 

Liley, D. T. J., Foster, B. L., and 
Bojak, I. (2012). "Cooperative 
populations of neurons: mean 
field models of mesoscopic brain 
activity," in Computational Systems 
Neurobiology, ed N. Le Novere 
(Netherlands: Springer), 315-362. 

Liston, A. D., Lund, T. E., Salek- 
Haddadi, A., Hamandi, K., Friston, 
K. J., and Lemieux, L. (2006). 
Modelling cardiac signal as a 
confound in EEG-fMRI and its 
application in focal epilepsy study. 
Neuroimage 30, 827-834. 

Logothetis, N. K., and Wandell, B. 
A. (2004). Interpreting the BOLD 
signal. Annu. Rev. Physiol. 66, 
735-769. 

Lopes da Silva, F. H., Hoeks, A., Smits, 
H., and Zetterberg, L. H. (1973). 
Model of brain rhythmic activity: 
the alpha-rhythm of the thalamus. 
Kybernetik 15, 27-37. 

Marreiros, A. C, Kiebel, S. J., 
Daunizeau, J., Harrison, L. 
M., and Friston, K. J. (2008). 
Population dynamics under the 
Laplace assumption. Neuroimage 
44, 701-714. 

Moran, R. J., Jung, E, Kumagai, T., 
Endepols, H., Graf, R., Dolan, 
R. J., et al. (2011). Dynamic 
causal models and physiological 
inference: a validation study 
using isoflurane anaesthesia in 
rodents. PLoS ONE 6:e22790. doi: 
10.1371/journal.pone.0022790 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



January 2013 | Volume 6 | Article 103 | 15 



Daunizeau et al. 



Predictive validity of sDCM 



Nunez, P. L. (1974). The brain wave 
equation: a model for the EEC 
Math. Biosci. 21,279-297. 

Penny, W., Joao, M., Flandin, 
G., Daunizeau, J., Stephan, 
K. E., Friston, K. J., et al. 
(2010). Comparing families of 
dynamic causal models. PLoS 
Comp. Biol. 6:el000709. doi: 
10.1371/journal.pcbi.l000709 

Pinotsis, D. A., Moran, R., and Friston, 
K. J. (2012). Dynamic causal mod- 
elling with neural fields. Neuroimage 
59, 1261-1274. 

Roebroeck, A., Formisano, E., and 
Goebel, R. (2011). The identifi- 
cation of interacting networks in 
the brain using fMRI: model selec- 
tion, causality and deconvolution. 
Neuroimage 58, 296-302. 

Rosa, M. J., Daunizeau, J., and Friston, 
K. J. (2010a). EEG-fMRI integra- 
tion: a critical review of biophys- 
ical modelling and data analysis 
approaches. /. Integr. Neurosci. 9, 
453-476. 

Rosa, M. J., Kilner, J., Blankenburg, 
R, Josephs, O., and Penny, W. 
(2010b). Estimating the transfer 
function from neuronal activity to 
BOLD using simultaneous EEG- 
fMRI. Neuroimage 49, 1496-1509. 

Schmuel, A. (2010). "Locally mea- 
sured neuronal correlates of 
functional MRI signals," in EEG- 
fMRI- Physiology, Technique 
and Applications, eds C. Mulert 



and L. Lemieux (Berlin, 
Heidelberg: Springer), 63-82. 
ISBN: 978-3-540-87918-3 

Schultze-Kraft, M., Becker, R., 
Breakspear, M., and Ritter, P. 
(2011). Exploiting the potential of 
three dimensional spatial wavelet 
analysis to explore nesting of 
temporal oscillations and spatial 
variance in simultaneous EEG- 
fMRI data. Progr. Biophys. Mol. Biol. 
105, 67-79. 

Seung, H. S., Richardson, T. J., 
Lagarias, ]. C, and Hopfield, J. J. 
(1995). "Minimax and Hamiltonian 
dynamics of excitatory- inhibitory 
networks," in NIPS y 97 Proceedings 
of the 1997 Conference on Advances 
in Neural Information Processing 
Systems, Vol. 10 (Cambridge, MA: 
MIT Press), 329-335. 

Smith, S. M., Miller, K. L., Salimi- 
Khorshidi, G., Webster, M., 
Beckmann, C. R, Nichols, T. 
E., et al. (2011). Network modelling 
methods for fMRI. Neuroimage 54, 
875-891. 

Sporns, O. (2010). Networks of the 
Brain. Cambridge, MA: MIT Press. 
ISBN: 9780262014694. 

Srinivasan, R., Tucker, D. M., and 
Murias, M. (1998). Estimating 
the spatial Nyquist of the human 
EEG. Behav. Res. Methods lustrum. 
Comput. 30, 8-19. 

Stephan, K. E., Kasper, L., Harrison, 
L., Daunizeau, J., Den Ouden, 



H., Breakspear, M., et al. (2008). 
Nonlinear dynamic causal mod- 
els for fMRI. Neuroimage 42, 
649-662. 

Stephan, K. E., Penny, W. D., 
Daunizeau, J., Moran, R. J., 
and Friston, K. J. (2009). Bayesian 
model selection for group studies. 
Neuroimage 46, 1004-1017. 

Stephan, K. E., Weiskopf, N., Drysdale, 
P. M., Robinson, P. A., and Friston, 
K. J. (2007). Comparing hemo- 
dynamic models with DCM. 
Neuroimage 38, 387-401. 

Tononi, G., Sporns, O., and Edelman, 
G. M. (1994). A measure for 
brain complexity: relating func- 
tional segregation and integration 
in the nervous system. Proc. 
Natl. Acad. Sci. U.S.A. 91, 
5033-5037. 

Valdes-Sosa, P., Roebroeck, A., 
Daunizeau, J., and Friston, K. 
(20 11). Effective connectivity: 
influence, causality and biophys- 
ical modelling. Neuroimage 58, 
339-361. 

Vaudano, A. E., Laufs, H., Kiebel, S. 
J., Carmichael, D. W., Hamandi, K., 
Guye, M., et al. (2009). Causal hier- 
archy within the thalamo-cortical 
network in spike and wave dis- 
charges. PLoS ONE 4:e6475. doi: 
10.1371/journal.pone.0006475 

Von Holdt, R. E. (1965). Rational 
power of power series. Am. Math. 
Mon. 72, 740-743. 



Wilson, H. R., and Cowan, J. D. (1972). 
Excitatory and inhibitory inter- 
actions in localized populations 
of model neuron. Biophys. J. 12, 
1-24. 

Wilson, H. R., and Cowan, J. D. (1973). 
A mathematical theory of the 
functional dynamics of cortical and 
thalamic nervous tissue. Kybernetik 
13, 55-80. 

Conflict of Interest Statement: The 

authors declare that the research 
was conducted in the absence of any 
commercial or financial relationships 
that could be construed as a potential 
conflict of interest. 

Received: 26 March 2012; accepted: 31 
December 2012; published online: 18 
lanuary 2013. 

Citation: Daunizeau }, Lemieux L, 
Vaudano AE, Friston KJ and Stephan KE 
(2013) An electrophysiological valida- 
tion of stochastic DCM for fMRI. Front. 
Comput. Neurosci. 6:103. doi: 10.3389/ 
fncom.2012.00103 

Copyright © 2013 Daunizeau, 
Lemieux, Vaudano, Friston and 
Stephan. This is an open-access article 
distributed under the terms of the 
Creative Commons Attribution License, 
which permits use, distribution and 
reproduction in other forums, provided 
the original authors and source are cred- 
ited and subject to any copyright notices 
concerning any third-party graphics etc. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



January 2013 | Volume 6 | Article 103 | 16 



Daunizeau et al. 



Predictive validity of sDCM 



APPENDIX 1: APPROXIMATING THE LOCAL NEURAL FIELDS WITH A SPATIAL EIGENMODES EXPANSION 

Let V 2 = d 2 /dx 2 + 9 2 /9y 2 be the 2D-Laplacian operator defined on a bounded (square) 2D domain of size ||D||, where {x, y] = r 
are the two cartesian coordinates on the domain. Recall that eigenvectors w^' () (r) of V 2 are such that: 



= w (k) (x)w (1) (y) 



(Al) 



where we have used index pairs (k, I) for spatial frequencies along the x- and y-axis, respectively. The associated eigenvalues are given 
by: Xjt j = — (k 2 + I 2 ) Tt 2 /||D|| 2 . By convention, w^ 1 refer to the eigenfunctions of the lD-Laplacian, and its complex conjugate will 
be denoted as w^ k \ Now, projecting Equation (5) onto the (complex conjugate) spatial eigenmodes yields the eigendynamics of the 
neural field (Daunizeau et al., 2009a): 
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where the projections q^, ;(0 of the firing rate form a complete basis function set. 

Together with a truncated eigenmode expansion, Equation (A2) can be used to approximate the solution to Equation (5) of the 
main text. Such a truncation basically acts as a low-pass filter, neglecting the transient effects of higher spatial frequency modes. This 
can be motivated by noting that ( 1 ) the eigenmodes w' ' (r) of the Laplacian operator are of increasing spatial frequency as the order 
(fc, I) increases and (2) the associated eigenvalues are negative (except for Xoo = 0) with increasing magnitude. This implies that 
the dynamics of high-order spatial modes will be strongly damped and thus quickly disappear in the absence of external forcing. 

We now focus on deriving an approximation to the projections 5^ /(f) of the firing rate. 

First, let us Taylor-expand the sigmoidal mapping of PSP activity around zero: 



S(T|(r, 0) 
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m = 0 



where = |^| is the /nth derivative of the sigmoidal mapping evaluated at r| = 0 and r\"'(r, t) is the mth power of mean 

PSP r|(r, f). Assuming that Equation (Al) holds for the eigenvectors of the Laplacian operator defined onto the cortical domain D, 
Equation (6) becomes: 



Ak) 



k I 

= J2~ z *\t,y)[ W M(x)] k 

k 



(A4) 



where z\ (t, y) is a dummy eigenmode dynamics, that is still a function of space (more precisely, of the second coordinate/). Inserting 



Equation (A4) into (A3) yields the following power series of power series in w 



in. 



00 / ,\> 

m = 0 m ' \ k / 

00 

= T, ^S™J2A k (t,y,m)wM(x) 



(A5) 
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where the coefficients A k (t, y, m) obey the following recursive equation (Von Holdt, 1965): 



A k (t,y, m) = -i- Td =1 ((m+ l)i - k)zf {t)A k -i{t,y, m) 

r % r (A6) 



The contribution Ti-(f , y) of the fcth spatial mode to the expansion in Equation (A5) is given by: 

S(T)(r, 0) = ^T t (t,y)w w (x) 



(A7) 

r k (t,y)= J2 —s 0 lm] A k (t,y,m) 



These can be analytically derived from the recursive relationship given in Equation (A6). They are given in Table Al bellow, up to 
order k = 3. 

It can be seen form Table Al that, except for the fundamental mode, the spatial mode contribution T k is of the form: 



T,«^o([5 1 y)sw(5S 0) ) 



;'=i 
k 



= ^o([i 1 F)p^(i< 0) ) n (l-ms(sj°>)), (A8) 

j=l m=l 

where O (x 1 ) is a polynomial of order j in x and p is the sigmoidal slope. The second line in Equation (A8) follows from the fol- 
lowing property of the sigmoid function: S^\x) = p'S(x) Yl m = i0- ~ rnS(x)) and the last line would follow from a first-order Taylor 
expansion of the sigmoidal mapping. This means that the lower order term in each spatial mode contribution T k (t,y) is already a 
polynomial of order three in the dummy eigenmodes dynamics z\ ■* (t, y). We will therefore truncate the contribution to its first term, 
i.e.,: Yrj(f, y) = S (zf\t, y)^j for k = 0, and T k (t, y) ~ zf\t, y)S' (zf\t, y)j for k > 0. This yields the following approximation to 
the firing rate field: 

S(t!(r, t)) ~ s(zf\t,y)) w (0) (x) + £ P S' (2< 0) (f,y)) ifV,y)w (t) W (A9) 

it>0 



Table A1 | This table gives both the coefficients A k (t. y, m) (left) and the contribution of the ftth spatial mode to the expansion in Equation 
(A4), as given by T k (t, y) (right) up to k = 3. 

A k [t, y, m) T k (t y) 



k = 0 



k = 2 



* = s m - z f > [zf ] m - + aw + w h <1) ] [ 2 T ^ (3>s ' ( J ! 0) ) + W 5 " (^ <0) ) + i- m sl3] 
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Equation (A9) can now be used to approximate the integral with respect to x in Equation (A2): 



dy 



Sk,0) = / ^~ n (y) f w ( - k) {x)S^(r,t))dx 

~hf w { - l) (y)s{zf ) (t,y)) dy+(l- S k )p J w ( - ,) (y)S'(zf\t,y))z (k \t,y)dy (A 10) 

= hj w^(y)s(zf\t,y))dy+(l-8 k )pJ24 J (t) f w (l '-'\y)S' '(if 1 \t,y))dy 

where 8^ is a Dirac delta, which is zero except when k = 0. 

Similarly, the sigmoidal terms in the integrands of Equation (A10) can be approximated using power series of power series, as 
follows: 



S (zf\t,y)) « S (zf' ■»(*)) + £ pS' (zf' °\t)) zf' \t)^\y) 

S' (zf\t,y)) « S' (z< 0 ' 0) ( f) ) w(°)(y) + £ pS" (z?' 0) ( t) ) zf' \f)«fl(y) 



(All) 



;>o 



Again, inserting Equation (All) into (A10), and ignoring high-order terms 1 yields: 

q k j(t) = hhS (zf' 0) (f)) + (1 - hh) PS' (zf' 0) (f)) zf- "(f) 



Inserting Equation (A12) in (A2) yields: 



zf'V) = 



±(G°s(zf- 0 \t))-4 0 -°\t)) iffc=i=0 
£ ((a 2 X w - 1) zf °(r) + P G°S' (zf ' °V)) zf' ;) (t)) otherwise 



(A12) 



(A13) 



The second line of Equation (A13) is Equation (8) of the main text, having re-ordered the eigenmodes index pairs (k, I) — >■ /c, for the 
sake of notational simplicity. 

The critical properties of this derivation can be summarized as follows: 

First, the fundamental eigenmode response will be slower than higher-order eigenmodes (because of smaller dampening terms). 
Second, the response of fast eigenmodes is enslaved by the slow eigenmode. Third, the fast eigenmodes do not feedback onto the slow 
eigenmodes. Taken together, we can think of the neural field as a feedforward system: inputs to the field activate the slow eigenmode 
zf' °^ (f), which controls the fast modes z\ ' (> (f)- This is further discussed in the main text of this document. 

APPENDIX 2: DERIVING THE FREQUENCY SPECTRUM OF EIGENMODES 

Recall that the eigenmode dynamics obey the following ODE (cf. Equation 8 in the main text): 

i(k) _ *_(*:) 



Ax w + v 



:(zf) -£(l-o»X*) 0 



_ 2 

T —I 



(A14) 



where v = 0(z 2 ) are high-order error terms that act as a perturbation to the deterministic evolution of eigenmodes. As a simplifying 
assumption, we consider that v behaves as i.i.d. white noise, the Fourier transform F [z®] of Equation (A8) writes: 



F Tz w j = (2nrwl 3 - A)" 1 F[v] 

OC (HtZwIt, — A) -1 I3 

where 1 3 is a 3 x 1 vector of ones, and I3 is the 3x3 identity matrix. 



(A15) 



1 These are O (z 2 ) terms that derive from the second term in the right-hand side of the second line of Equation (All) 
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Note that the inverse translated Jacobian writes: 



(2z'jtwI3 — A) 



1 



|2nta)l3 — A| 



l=k(io) (fto + \) 



!*(«>) 



c(z< 0) )(ico+§) ja,(io+2) + ^ c( z <°>) 



(A16) 



£*(©) = «»+ — (l-a 2 X t ) 
2a 

where ^jt(a)) is the only term that depends upon the mode order k, and the determinant of the translated Jacobian is given by: 



|2ijtcoI 3 - A| = ito^fc(oo) h'co + — (zS 0) ) + %k(<*>)-^ 



= %k(w) ico ico + - + 



(A17) 



where ^(w) is given in Equation (A12) (with appropriate eigenmodes index reordering (fc, Z) — > k). The explicit Fourier transform of 
eigenmode dynamics now derives from inserting Equation (A16) into (A15): 



1 



|2mool3 — A| 



% k (w) (tea + | + 1) + * 
(i» + |)( C (^>Ui») + ^ +C (^) 
^( W )(i W -i) + ^(; W + £ :(z< 0) )) 



(A18) 



Now, we are interested in the total frequency power of EEG signal sampled at some position r' on the scalp that is induced by mean 
membrane depolarization: 



P(w, r') 



f L(r, ^FMt, t)]dt 
Jd 

k 

a (k) (r') = I L(r,r')w ik \r)dr 
Jd 



(A19) 



where the second line derives from the eigendecomposition of the neural field (cf. Equation 6), L (r, r') is a spatial convolution operator 
that describes the quasi-static volume propagation of electromagnetic fields through head tissues, and a® (r') captures the contribu- 
tion of the fcth eigenmode at position r' on the scalp. It is known that the geometrical properties of electrical conductivity within head 
tissues eventually shape the EEG susceptibility to different cortical spatial scales. For the sake of simplicity, we will assume the skull 
acts as a perfect low-pass filter (see e.g., Srinivasan et al., 1998), which induces a cut-off order k c , i.e.,: 



a ( V) = { 



a (0 V) itk<k c 
0 if k > k c 



(A20) 



Equation (A20) basically neglects the continuous loss of sensitivity as one increases the spatial frequency of (cortical) membrane 
depolarization fields. Inserting Equation (A20) into (A19) yields Equation (9) of the main text, i.e.,: 



P(oo,r') = a (0) (r') 2 



k<kc 



„(0) 



(rV £ 



% k (a>) (ico+ 2 - + l) + 



k<k c 



a^)(-(- + f) + i 2 )-" C H 0) ) 

where a' 0 ' (r') is a scaling factor that will be different for different EEG sensors. 



(A21) 
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Equation (A21) is a valid approximation of the EEG power spectrum at position r' on the scalp generated by one local brain region. 
The generalization to a network of remote brain regions requires the specification of the relative contribution of each region, which 
depends on their relative distance with EEG scalp sensors. In principle, this can be derived from knowledge of the forward operator 
L (r, r') (see e.g., De Munck, 1988). This is beyond the scope of the present study. 

Lastly, note that the numerical values of the Laplacian's eigenvalues are approximated from the 2D-Euclidean case, as follows: 
Xjt, i = — (k 2 + I 2 ) Jt 2 /||D|| 2 . The effective frequency cut-off k c is such that k 2 + I 2 < k 2 . In practice, we use k 2 = 16 and ||D|| = 1 cm. 
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